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Abstract 

While it is widely believed that gravity should ultimately be treated as a 
quantum theory, there remains a possibility that general relativity should 
not be quantized. If this is the case, the coupling of classical gravity to 
the expectation value of the quantum stress-energy tensor will naturally 
lead to nonlinearities in the Schrodinger equation. By numerically in- 
vestigating time evolution in the nonrelativistic "Schrodinger-Newton" 
approximation, we show that such nonlinearities may be observable in 
the next generation of molecular interferometry experiments. 
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1. Introduction 

The first attempts to quantize general relativity appeared in the early 1930s, and in 
the seventy years that have followed, we have learned a great deal about quantum field 
theory and gravity But despite the hard work of many outstanding physicists, a complete, 
consistent theory of quantum gravity still seems distant [1]. Given the severe difficulties, 
one might reasonably ask whether the whole project could be a blind alley. Perhaps our 
prejudice that everything in Nature should be quantized is simply wrong; perhaps general 
relativity, a theory of spacetime, is fundamentally different from theories of fields within 
spacetime. 

Somewhat surprisingly, the possibility that gravity is essentially classical has not yet 
been excluded. The simplest model of classical general relativity coupled to quantum 
matter, sometimes called "semiclassical gravity," was proposed forty years ago by M0ller [2] 
and Rosenfeld [3]. The Einstein field equations become 



where the operator-valued stress-energy tensor is replaced by an expectation value. As a 
Hartree-like approximation to quantized gravity, such a system certainly makes sense. But 
as Kibble and Randjbar-Daemi stressed [4], viewed as a fundamental theory, such a model 
leads to nonlinearities in quantum mechanics: the Schrodinger equation for the wave func- 
tion \tp) depends on the metric, which in turn depends, through (jl.lj) . on While more 
complicated models are possible, any such theory must couple classical gravity to quantum 
sources, and it is hard to see how to do so without introducing similar nonlinearities. 

While the literature contains a number of criticisms of semiclassical gravity [4,6-9], none 
seems decisive [3,5,10,11]. For example, one might argue that measurements with nonquan- 
tized gravitational waves could violate the uncertainty principle for quantized matter [7]. 
But there are intrinsic limitations to the measurement of even a classical gravitational 
field [5, 12]; it is possible, for instance, that the necessary measurement would require 
an apparatus massive enough to collapse into a black hole [11]. On the experimental side, 
neutron interferometry [13] and microscopic deflection experiments [14] show that quantum 
matter interacts gravitationally as expected, but these results do not require quantization 
of the gravitational field itself. More direct experimental tests have been suggested using 
superpositions in Bose-Einstein condensates [15] or, in principle, gravitational radiation 
from quantum systems [16], but these are not yet practical. The standard mechanism for 
the appearance of inhomogeneities in inflationary cosmology is probably incompatible with 
semiclassical gravity, requiring at least linear quantum fluctuations in the metric [17], but 
there are certainly other ways to produce the observed inhomogeneities. 

Semiclassical gravity is probably excluded in a strict Everett interpretation of quantum 
mechanics [8], but one would like to have a result that does not depend on measurement 

*Dirac was also aware of this: see [5], p. 1. 
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theory or the interpretation of quantum mechanics. It is therefore important to ask whether 
the nonlinearities in semiclassical gravity might be experimentally accessible. 

The full semiclassical equations (jl.lj) are extremely difficult to analyze. But for many 
purposes, the nonrelativistic Newtonian approximation should be sufficient. We are thus 
led to the ''Schrodinger-Newton equation" [18,19], 

ih?¥- = — — V 2 </> - m<S>ip, V 2 $ = AnGmM 2 , (1.2) 
at 2m 

a Schrodinger equation for matter coupled to a classical gravitational potential that has as 
its source the expectation value of the mass density. This system has been studied in the 
past [20-23], and a fair amount is known about the lowest eigenf unctions and eigenvalues, 
but time evolution remains much more poorly understood [24-26] . 

In the remainder of this paper, we will investigate the Schrodinger- Newton evolution of 
an initial Gaussian wave packet. Qualitatively, we will show that self-gravitation can slow 
or stop the spreading of the wave packet. We find that despite the weakness of gravity, the 
resulting suppression of interference may be observable in the next generation of matter 
interferometry experiments [27]. 



2. Setting up the problem 



Let us first consider a few analytic properties of the Schrodinger-Newton equation. 
Note that although the equation is nonlinear, time evolution preserves the norm of if), and 
a conserved probability current can be written down: 



^1 
dt l 



ih 
2m 



(2.1) 



Thus, as in standard quantum mechanics, we can interpret \ip\ 2 as a probability density. 
We will be interested in an initial Gaussian wave function 



V(r,0) 



a\3/4 
n) ' 
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(2.2) 



with width a 1//2 . In principle, we expect a two-parameter family of solutions, labeled by 
a and m. But the Schrodinger-Newton equation is invariant under the rescaling 



m 



jjm, 



x 



(2.3) 



so if ip(a,m; x,t) is a solution, so is /i 9//2 ?/>(/i 6 a, fim; /i~ 3 x, /i 5 t). It therefore suffices to 
consider a one-parameter family of solutions. 

We have investigated the dynamical Schrodinger-Newton equation (jl.2j) with initial 
condition ()2.2j) numerically, using a homegrown PDE solver specifically designed for this 
problem, along with associated tools to analyze and organize the data. Details may be 
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found in P. J. Salzman's dissertation [28], and will appear in a subsequent paper. Briefly, 
the Schrodinger equation can be discretized by considering the action of a time evolution 
operator on the known wave function at some time t to yield the unknown wave function 
at some later time t' . For example, two common schemes are the "explicit method," which 
translates the wave function forward in time, 

ip(r,At) =e- iitAt/h i)(r,0), 

and the "implicit method," which translates the wave function backward in time, 

e +iAAt/h ^{r,At) = V>(r,0). 

Here H is the Schrodinger-Newton Hamiltonian, which involves both derivatives and the 
integral 

all space 

the solution to the Poisson equation (jl.2j) . With either scheme, one can perform a Taylor 
expansion of the energy operators to any desired order, requiring successively higher order 
derivatives that can be approximated and converted into a numerical algorithm. 

Each of these schemes, however, has properties which make it undesirable to use for 
a numerical PDE solver. The explicit scheme, while simple and calculationally inexpen- 
sive, can be shown to be numerically unstable. The implicit scheme, while numerically 
stable, requires that we find an inverse operator, which is complicated and calculationally 
expensive. Additionally, neither scheme is unitary. 

We chose to use Cayley's form [29], which is an average of the explicit and implicit 
methods: 

e iHAt/2h^ At / 2 ) = e -iHAt/2H^ _ At / 2 ). (2.5) 

Cayley's form is unitary and is of higher order than either the implicit or explicit method. 
The evolution operator on either side of Cayley's form was Taylor expanded to first order, 
giving a second order accurate algorithm. Derivative terms were numerically approximated, 
and yielded a tridiagonal system of linear equations from which the wavefunction at the 
next timestep could be computed. The integral I$n technically involves the wave function 
at the "current timestep," at which it is unknown; we linearized the integral by using 
the known wave function at the previous timestep. Accuracy was ensured by successively 
decreasing the timestep At and looking at the limiting behavior as At approached machine 
precision. 

A wide assortment of techniques was employed to further test the accuracy of our 
method. As one important benchmark, for each set of input parameters the program was 
rerun with the potential turned off, and the results were compared with the exact analytic 
solution for the free particle; the differences were found to be negligible. 
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Figure 1: Comparison of Schrodinger-Newton and free particle wave functions. For intermediate 
masses (a), the evolution becomes complex; for larger values (b), the wave function "collapses." 



3. Numerical Results 

The Schrodinger-Newton equation was repeatedly solved using a fixed wave packet 
width a = 5 x 10 16 m~ 2 while varying the mass m. The results may be summarized as 
follows (with masses in unified atomic mass units): 

1. For small masses (< 4.5 u), the behavior is essentially indistinguishable from that of 
a free particle. As m increases, the wave packet spreads more slows, as one might 
expect from "self-gravitation." 

2. For masses between 1.9 x 10 3 u and 7.2 x 10 3 u, the behavior is complex; the wave 
packet typically fluctuates rapidly and develops growing oscillations (figure la). This 
seems to be related to the behavior found in [24-26]. 

3. For masses between 7.8 x 10 3 u and 2.9 x 10 14 u, the wave packet "collapses," shrinking 
in width (figure lb). 

4. For larger masses, the wave packet appears stationary; we have not been able to run 
the program long enough to determine the behavior. 

One can obtain a crude estimate of the "collapse" mass by noting that for a free particle, 

/ 2 2 \ 1/2 

the peak probability density of a Gaussian wave packet occurs at r p ~ a' 1 ' 2 ( 1 + ^rt 2 J , 

"accelerating" at a speed r p ~ h 2 /m 2 r p 3 . Equating this with the acceleration due to gravity 
at t = yields a mass 

m ~ (^r) V3 ~ 1010 u ' (3 - 1} 
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Figure 2: "Phase diagram" for Schrodinger-Newton solutions. Region A: wave packets spread; 
region B: complex behavior; region C: collapse; region D: undetermined by our simulation. 

lying roughly at the middle of the mass range for which we see "collapse." This expression 
behaves properly under the scaling (|2.3|) . and could be guessed by dimensional analysis; a 
key result of our numerical simulations is that onset of "collapse" occurs at significantly 
smaller masses, presumably reflecting the nonlinearity of the evolution. 

Using the scaling ()2.3|) . we can summarize our results by the "phase diagram" shown in 
figure 4. Of particular interest for experiment is the "collapsing" phase (C). As the figure 
illustrates, the wave packet of fixed width w = a _1//2 will shrink in width if its initial mass 
lies within a range m_(w) < m < m + (w). Numerically, with w in microns and m in unified 
atomic mass units, we find 

m_/l u = 1300(w/l fim)- 1/3 , m+/l u = 4.8 x 10 13 (w/l /WT 1/3 . (3.2) 

The characteristic collapse times in nanoseconds, obtained numerically and scaled according 
to (E3D, are 

T_/lns = 1.2 x 10- 4 (w/lfim)- 5/3 T+/1 ns = 1.2 x Kr 2 (w/1 /im)" 5/3 (3.3) 

4. Experimental tests 

A "collapsing" wave packet will lead to a suppression of interference, which should in 
principle be observable. Most recent experiments in matter-wave diffraction have used 
Talbot-Lau interferometry [30], in which an image of a diffraction grating with slit spacing 
d appears at a distance = d 2 /X from the grating. The heaviest molecule for which 
interference has been observed to date is fluorofullerene, C 6 oF 4 8, with a mass of 1632 u [31]. 
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The grating slits in this experiment have a width w ~ .5 /mi. From (J3.2)) . semiclassical 
gravity would predict a loss of interference for a wave packet of this width for masses 
greater than m ~ 1600 u; C60F48 lies just at the edge of this range. 

This is, of course, an oversimplification: the molecular wave packets in [31] are not 
spherically symmetric Gaussians. The grating slits limit the width in one direction, say y, 
while the widths in the x and z directions are determined by other factors. Repeating the 
argument that led to (|3.1j) for a cylinder and a slab, we might expect a new limiting mass 
on the order of 

m'_ = m_(w x /w) 1/3 (w z /w) 1/3 (4.1) 

In the fluorofullerene experiment, w z was controlled by a height limiter with a width of 
100 fim, giving a factor of about 6 in (|4.1|) . The appropriate value for the width in the 
direction of the beam is less clear. As a pessimistic estimate, we note that the distance 
from the first grating, which is responsible for the transverse coherence of the beam, and 
the second grating, responsible for the interference, was .38 m; using this value for w x in 
(|4.1|) would give a factor of about 90. Thus while the results for fluorofullerene interference 
probably do not directly probe semiclassical gravity, they appear to come within two to 
three orders of magnitude of a real test. 

To proceed further, we need both experimental and theoretical work. On the theory side, 
simulations should be done with more realistic wave packet profiles. It is also important to 
see how sensitive our results are to the exact (Gaussian) form of the initial wave function. 
In particular, since a stable spherically symmetric ground state exists [24], some initial 
profiles must have different behavior than what we have seen. 

On the experimental side, some progress can come from using narrower slits and from 
controlling the longitudinal width w x , perhaps with a shutter to restrict the time each 
molecule enters the apparatus. The most useful gain, though, would probably come from 
measurements of higher mass molecules. This is not an unreasonable goal: experimentalists 
have argued that by using optical gratings, it may be possible to exhibit interference for 
molecules with masses of up to 10 6 u [27,30,32,33]. If the next generation of matter wave 
interferometers can come even close to this limit, an unambiguous test of semiclassical 
gravity should be possible. 
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